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ABSTRACT 

A model is developed for the transport of heat and solute in a system of double-diffusive layers under astrophysical 
conditions (viscosity and solute diffusivity low compared with the thermal diffusivity) . The process of formation of the 
layers is not part of the model but, as observed in geophysical and laboratory settings, is assumed to be fast compared 
to the life time of the semiconvective zone. The thickness of the layers is a free parameter of the model. When the energy 
flux of the star is specified, the effective semiconvective diffusivities are only weakly dependent on this parameter. An 
estimate is given of the evolution of layer thickness with time in a semiconvective zone. The model predicts that the 
density ratio has a maximum for which a stationary layered state can exist, R p ^ Le -1 / 2 . Comparison of the model 
predictions with a grid of numerical simulations is presented in a companion paper. 
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1. Introduction 

In stellar evolution, 'semiconvection' denotes the situa- 
tion where a thermally unstable stratification is stabilized 
against (adiabatic) overturning by a gradient in composi- 
tion (called solute in the following; typically the Helium 
concentration, increasing with depth). It was first rec- 
ognized as a complication in the calculation of stellar 
structure by R. J. Tayler (1953). Uncertainty whether 
the stabilizing gradient should be ignored ('according to 
Schwarzschild'), included in the condition for overturning 
('according to Ledoux'), or something in between has led 
to a number of different recipes for mixing of the Helium 
gradient. The evolution of the star subsequent to the semi- 
convective phase is sensitive to these differences (e.g. Weiss 
1989, Langer et al. 1989, Langer 1991, Stothers & Chin 
1994, Langer and Maeder 1995). The fluid mechanics en- 
countered in geophysics in the same case of a thermally 
unstable stratification stabilized by a stabilizing solute is 
called double diffusive or thermohaline convection. 

The observations show that such a stratification forms a 
stack of many thin layers, called a 'staircase', each consist- 
ing of overturning convection sandwiched between stable 
steps in composition and temperature (Turner & Stommel 
1964, Padman & Dillon 1987, Schmid et al. 2010 and ref- 
erences therein). In effect, this is a 'compromise between 
Schwarzschild and Ledoux'. Correspondingly, the net trans- 
port coefficients (of heat and solute) are intermediate be- 
tween those of convection and diffusion. 

The reason for this layer formation are understood (for 
references, cf. Spruit 1992, hereafter S92, and Zaussinger & 
Spruit 2013, hereafter ZS13). The main theoretical contri- 
bution in this context has been the work of Proctor (1981), 
who showed analytically that a finite amplitude layered 
state exists for conditions when the stratification is still 
linearly stable: layering is the result of a subcritical insta- 
bility. His analysis applies to the case of vanishing diffusiv- 
ity of the solute, and Prandtl number not exceeding 0(1), 



which is the astrophysically relevant case. In this limit, lay- 
ered states exist whenever the Rayleigh number exceeds 
the critical value for normal convection, independent of the 
strength of the stabilizing solute gradient. An energy argu- 
ment (S92) shows that the layered state in this case can be 
reached with an initial perturbation of vanishing amplitude 
as the layer thickness decreases. 

Attempts have been made to address the astrophysi- 
cal problem with direct numerical simulations (Merryfield 
1995, Biello 2001, Rosenblum et al. 2011). This encounters 
two problems: the very high thermal diffusivity in a stel- 
lar interior (very small Prandtl number) cannot be matched 
without some form of approximation. More importantly the 
quantity of main interest, the effective mixing rate, depends 
on the thickness of the double diffusive layers formed, a 
quantity that is not a stable outcome of the simulations. 

It makes sense to disentangle the 'semiconvection prob- 
lem' into two parts: on the one hand the physics that de- 
termines the thickness of the layers, on the other hand the 
effective mixing rate for a given layering state. This sepa- 
ration is especially meaningful because observations in geo- 
physics and laboratory experiments show layer thickness 
to be a slowly changing function of time compared with 
the overturning times within the layers. In the following we 
concentrate on the second question, that is, the mixing rate 
is studied as a function layer thickness. 

It turns out that the simple observation of a layer struc- 
ture consisting of an overturning zone between stable zones 
is sufficient to derive a predictive model for the effective 
transport coefficients (Sect. [3|. It is sufficiently quantita- 
tive to be tested against the results from numerical simu- 
lations. The translation of the model to the astrophysical 
case of a compressible stratification can be done exactly in 
the limit where the layers are thin compared with the pres- 
sure scale height. It is given in sec. [6j and yields an easily 
implcmcntable prescription for the mixing rate. 

In the present treatment, the transition in composition 
and temperature between neighboring layers has a finite 
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Fig. 1. Notation used, showing temperature as a function of 
depth through half of a semiconvective layer of thickness d, 
from the middle (z = 0) of the stagnant zone (of thickness 
d s ), to the middle of the overturning zone at z = d/2. 
AT is the temperature difference across the whole layer, 
AT across the overturning zone. The solute profile has the 
same shape but different amplitudes AS and AS. 



width, as opposed to S92, ZS10, where it was approximated 
as a step. This generalization turns out to make no signif- 
icant difference for practical astrophysical application, but 
is essential for meaningful comparison of the theory with 
numerical results. Section [7] finally gives a (less quantita- 
tive) estimate of the layer thickness and its evolution in 
time. 

The astrophysical term 'semiconvection' will be used in 
the text interchangeably with the terms 'double diffusive' 
or 'thermohaline' convection used in laboratory and geo- 
physical literature^] 

2. Physics of a semiconvective layer 

2.1. Notation and definitions 

As in the above, we will use the term 'layer' for an in- 
dividual double-diffusive step in the semiconvective stair- 
case, and the term 'solute' for the stabilizing component 
(e.g. Helium in Hydrogen). Such a layer consists of a zone 
of overturning convection between adjacent stagnant zones. 
In the stagnant zones overturning convection is suppressed 
by a stable (TV 2 > 0) density gradient, and transport of 
heat and solute takes place by diffusion. 

As illustrated in Fig. [l] depth through the layer is 
counted from the middle of a stagnant zone. The thick- 
ness of the layer as a whole is d, that of the stagnant zone 
d s . The temperature difference across the layer is AT, the 
solute difference AS. The stabilizing influence of the so- 
lute is expressed conveniently in terms of the density ratio, 
the ratio of density differences caused by temperature and 
solute differences AT and AS across the layer: 



R p = /3 AS /a AT, 



(1) 



1 Existing nomenclature is somewhat ambiguous. The oppo- 
site case of a destabilizing composition gradient in a stable ther- 
mal gradient is called 'saltfingering', but in geophysics also 'ther- 
mohaline convection', where our semiconvective case is often 
called 'diffusive convection' (e.g. Schmitt 1994). 'Double diffu- 
sive' is usually meant to cover both cases. 



where a, the thermal expansion coefficient, is the relative 
density decrease per unit temperature and /3, called the ha- 
line contraction coefficient, is the relative density increase 
per unit of solute concentration. Since density increases 
with S and decreases with T, the density differences are 
of the opposite sign when AT and AS* both increase with 
depth (in the direction of gravity) as in our semiconvective 
case. 

The flow in the overturning zone of the layer is driven 
by a temperature difference AT (< AT) and opposed by a 
solute difference AS (< AS). Associated with these is an 
'internal density ratio', related to the overall density ratio 
by 

AS AT 



R,, 



pAS/aAT = R p — — . 
1 9 AS AT 



(2) 



Obviously, the values of AT and AS depend on where we 
define the boundaries between stagnant zone and overturn- 
ing zone. Since the internal layer is actually convecting, the 
boundary has to be set at a point in the T, S profile where 



R p < 1 (see also Sect.|3^ 

The overturning zone is characterized by a Rayleigh 
number; for convection in a layer of thickness D, with tem- 
peratures T t and Ti at top and bottom, the Rayleigh num- 
ber is 

Ba=g q(r b -r t )2 3 3 

where g is the acceleration of gravity, kt the thermal diffu- 
sivity, and v the (kinematic) viscosity. In our case (Fig. [T]) , 
D has the value d — d s . The critical value for onset of con- 
vection is of order Ra c = 1400 (for no-slip boundary condi- 
tions). If viscosity is low, (Prandtl number Pr= v/kt <C 1), 
the heat flux at large Rayleigh number becomes indepen- 
dent of viscosity, (a fact that is used implicitly in the 'mix- 
ing length' formalism for convection in a stellar interior). 
At high Ra, a more relevant quantity to characterize the 
heat flux in this case is the modified Rayleigh number Ra* : 



Ra* = Pr Ra = 



ga{T h - T t )D 3 



(4) 



Its square root can be read as the ratio of the thermal 
diffusion time D 2 / kt to the free fall time of the density 
contrast a(T> — T) over the distance D. 

Let F be the (time averaged) heat flux across the layer, 
and Td the flux in the absence of convection, i.e. when the 
temperature profile between T, and T t is determined by 
diffusion only. The Nusselt number is then defined as 

Nu T = F/F d . (5) 

Similarly, there is a Nusselt number for the solute flux F$: 

Nu s = T s /T sd . (6) 

In the absence of a solute, i.e. for normal (unstratified) 
laboratory convection, Nut is a function of Pr and Ra only 
(apart from boundary conditions and geometry). In the case 
of semiconvection, the Nusselt numbers are functions of two 
additional parameters characterizing the solute: the density 
ratio Rp and the Lewis number Le, the ratio of solute dif- 
fusivity ks and thermal diffusivity Kt- 



Lc = ks/kt. 



(7) 
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2.2. Transport coefficients: classical results 

The heat flux is determined by the boundary layers at 
top and bottom. As the overturning flow passes along the 
boundary, the temperature contrast with the boundary dif- 
fuses into the flow, and it is this temperature difference 
that carries the heat flux. If r is the time during which 
the flow is in contact with the boundary before descend- 
ing/ascending into the bulk, the depth Dt over which the 
temperature difference penetrates is Dt — (ktt) 1 / 2 . At 
high Rayleigh numbers, this depth is small compared with 
the layer thickness d. In a simple two-dimensional view of 
the flow the amount of fluid carrying this difference also 
scales with Dt- The same flow along the boundary de- 
termines how much solute contrast diffuses into the flow. 
The solute diffusion depth is thus D§ — (ks t ) 1 ^ 2 , an d the 
amount of solute flowing into the overturning zone is pro- 
portional to D$ and inversely proportional to the flow speed 
v along the boundary. 

The ratio of the convective fluxes of solute and heat is 
thus expected to scale as (ks/kt) 1 ^ 2 - On the other hand, 
the (diffusive) fluxes over the thickness of the overturning 
layer, in the absence of convection, scale as «s an d kt- 
In terms of the Nusselt numbers, this implies, in the limit 
Nut, Nus 3> 1 such that the boundary layers are thin com- 
pared with the layer thickness: 



Nu s = <jLe" 1/2 Nu T 



(8) 



where q is a numerical factor expected to be of order unity. 
If the heat flux is kept fixed, the effective diffusivity of the 
solute then scales as 



K s,cff 



= ftsNus = 9 (/cs^r?) Nut- 



(9) 



These are classical scalings, as found in previous analyses of 
double diffusive convection, e.g. in Turner (1980, 1985). In 
this form they do not address the dependence of the fluxes 
on the density ratio R p . The argument above implicitly as- 
sumes that all the density contrast in the solute boundary 
layer D$ is carried across with the flow. In reality, only a 
fraction of the boundary layer can flow across, namely the 
fraction that is not too buoyant to be carried down by the 
flow, respectively too heavy to be carried up. Taking this 
into account (Spruit 1992, hereafter S92), yields a correc- 
tion to j8j: 



Nu s = ^-Le~ 1/2 Nu T . (Nu T > 1) 
R„ 



(10) 



This scaling applies to the overturning zone; it would be 
valid for the layer as a whole only if the stagnant zone 
were absent. More precisely, it assumes that the thickness 
ds of the stagnant zone is smaller than the thickness of the 
boundary layers Dt,D$. This is not always the case (cf. 
Sect.|5]below). In the following, the presence of the stagnant 
zone is taken into account explicitly, by consistently taking 
the distinction between between Nusselt numbers for the 
overturning zone and those for the layer a s a whole into 
account. This results in a generalization of (10 1. 



2.2.1. 'Erosion' 



To determine the numerical factor q in ( 10 1 more quantita- 



tively than order unity, the hydrodynamic interaction of the 
flow with the stagnant zone has to be considered in more 



detail. This is somewhat beyond the scope of the model to 
be developed here, but we can identify a process involved, 
and use this to show that q is probably somewhat larger 
than unity. I will call q the 'erosion factor'. 

The unsteady flow in the overturning zone induces per- 
turbations in the stagnant zone. Since it is stably strati- 
fied, these take the form of internal gravity waves. These 
can contain internal structure on length scales (perpendicu- 
lar to the interface) less than the thickness of the stagnant 
zone. Diffusion of solute on these length scales increases 
the amount that has sufficiently low buoyancy to be car- 
ried with the overturning flow, hence we may expect q > 1. 
In the absence of a more detailed theory, its value can in 
principle be used as a fitting parameter, as long as it is not 
taken larger than order unity. In the following, however, I 
will ignore this option, and simply set 



q=l. 



(11) 



For astrophysical applications, the effective solute transport 
will turn out to be so low that tuning the erosion factor by 
order unity would have little effect anyway. At low Le, this 
erosion process takes place well inside the thermal bound- 
ary layer. It consequently affects only the transport of the 
solute; its effect on the transport of heat can be neglected. 

3. Model 

The layer of thickness d now consists (Fig. [T]) of a stagnant 
zone of thickness ds , and a overturning zone occupying the 
remainder of d. An incompressible (Boussinesq) fluid is as- 
sumed, to be generalized to the compressible astrophysical 
case in Sect. |6] 

In the stagnant zone, we assume that the transport of 
heat and solute is by diffusion only (i.e. ignoring the pos- 
sible 'erosion' effect discussed above). In the overturning 
zone of the layer, the flow is approximated as convection as 
it would take place in the absence of a solute. The presence 
of the solute affects the flow somewhat, but the amount of 
solute carried, and hence its influence on the flow, vanish 
in the limit of low solute diffusivity (S92, Schmitt 1994). 
Apart from the boundary conditions it sets, the stagnant 
zone has little effect on the flow inside the overturning zone. 
In the limit Le -C 1 considered here, the effect of the so- 
lute on the flow of heat in the overturning zone can thus be 
neglected. 

To describe transport of heat in the overturning zone, 
we use a fit from laboratory measurements for the Nusselt 
number as a function of the temperature difference. Since 
the thickness and temperature difference are different from 
those of the layer as a whole, the Rayleigh and Nusselt 
numbers of the overturning zone are distinguished here with 
a " . With the notation 



e = AT/ AT, 6 = d s /d, 



(12) 



they are related to Ra and Nut by (cf. eq. [4] and the geom- 
etry sketched in Fig. [I]) : 



Ra* = Ra* e(l - S) 3 , 
Nu T = Nu T (1 - S)/e. 



(13) 
(14) 



In the measurements of Niemela et al. (2000), at Ra up 
to 10 17 , the Nusselt number is well fit by a power law of 
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slope 0.309 and amplitude 0.124. The Prandtl number in 
these experiments is order unity, so Ra* and Ra are approx- 
imately the same. Combining the above, I approximate the 
heat flux in the overturning zone as given by 



Nun 



l + a(Ra* -Ra* c ) b , (a = 0.124, b = 0.309). (15) 



Since, as argued in 2.1 the dependenc e on Prandtl number 
is expected to be weak below Pr= 1, ( |T5"| will be assumed 
approximately valid for all Pr < 1 (tms can be checked 
with numerical simulations, see Zaussinger and Spruit 2013, 
hereafter ZS13). The constant 1 and the critical Rayleigh 
number have been added to better approximate the behav- 
ior at low Ra*. 

The buoyancy of the solute is opposite to that of the 
driving temperature difference, and the flow will only be 
driven by density differences of the unstable sign. This im- 
plies that the solute concentration contrast in the overturn- 
ing flow, S say, is limited by the temperature contrast T 
through 

RpS/AS < f/AT. (16) 

At the boundary between the overturning and the diffusive 
zones, these contrasts are S = AS, T = AT. Define then 
the location o f th is boundary as the point where the equals 
sign holds in (16 1, that is, the point where the buoyancy of 



the solute is just low enough for it to be carried with the 
thermally driven flow in the overturning zone (cf. discussion 
above). This yields the relation 



AS 
AS 



1 AT 
~Rp~AT' 



(17) 



In the overturning zone, the original scaling ([8| applies, 
i.e. Nug = Lc~ 1//2 Nut (assuming q = 1), since AS* is the 
density difference that can just be carried with the flow. 
At low Ra, Nug should approach unity at the same time as 
Nut- To accomplish this I adopt a slight modification: 



Nu s - 1 = IxT 1/2 (Nu T - 1 



(18) 



which is now assumed to hold for all Nut > 1- 

This completes the definition of the model. It defines 
the fluxes of heat and solute uniquely as functions of the 
external parameters. 

4. Effective diffusivities 

In a stationary state, the fluxes of heat and solute are con- 
stant with depth through the layer. The fluxes by diffusion 
in the stagnant zone are equal to the fluxes in the overturn- 
ing zone. Expressing this in terms of the Nusselt numbers, 
we first need to write Nut and Nus m terms of the Nusselt 
numbers for the layer as a whole. For temperature this is 
given by ( 14 ) : 

Nu T = Nut tz 6 ■ (19) 



(i-sy 

The equivalent relation for the solute is 



Nu s = Nu s 



(i-sy 



where (using 17) 



The fluxes in the stagnant zone are given by: 

Nu T - (1 - e)/S, Nu s = (1 - e s )/6. (22) 



(since carried by diffusion alone). With ( JT8| ), eqs. ( 19 )-( 22 1 
yield two equations for the unknown S and e: 

(l/5-l)(l/e-l)=Nu T , (23) 
(1/6 - l)(R p /e - 1) = 1 + IxT^Nut - 1), (24) 



with Nut given by (15). Eliminating Nut between these 
two yields 



1 - S - e/R p = (1-S- e)/Q, 



where 

Solving this for e: 



Q = i? p Le 1/2 . 



(l-S) 



1 - Q/R P 



(25) 
(26) 

(27) 



Since e must be a positive number, and R p is larger than 
1 (otherwise we would not be 'below Ledoux'), it follows 
that in a stationary state as envisaged Q must be less than 
1, or R p < Le -1 ' 2 . This is a necessary condition (inde- 
pendent of the Rayleigh number), but not sufficient. The 
actual, somewhat lower, valu e of the critical density ratio 
has to be determined from (fl9|)-(20) by solving for S as 
well as e. To see how this solution comes about, consider 
the left-hand and right-hand sides of ( 23 ) separately. They 



represent the heat flux in the stagnant and the overturning 
zones, respectively, and in a stationary state they are equal. 
Fig. [2] shows the tw o as functions of 6, with e taken from 
m\ and Nu T from (flil. 



The two intersection points of the curves are potential 
solutions for a steady state. To see which of the two is 
the relevant one, consider the slopes of the curves. At the 
equilibrium point with the smaller value of 6 the flux in 
the stagnant zone decreases more rapidly with S than the 
flux in the overturning zone. A small decrease of S away 
from the equilibrium would increase the flux in the stagnant 
zone relative to that in the overturning zone. This would 
result in a temperature deficit at the boundary between the 
two, causing the temperature to decrease there. This would 
reduce the heat flux in the stagnant zone again, the opposite 
of the assumed perturbation. This equilibrium point is thus 
the stable one; the intersect at the larger 5 is an unstable 
point. 

The necessary condition for existence of the layered 



state, R„ < Lc 



-1/2 



also figures prominently in Proctor 



e s = AS/ AS = e/R p 



(1981). In his analysis it shows up as a necessary condi- 
tion for its validity. 

4.1. Maximum density ratio 

The maximum density ratio can be determined as a func- 
tion of Ra* as the value for which the two curves in Fig. 
[2] just touch. The result is shown in Fig. [3] At large Ra*, 

(20) i? pmax (slowly) approaches the value Le" 1 ^. The behavior 
of Nut and S near R p max is shown in Fig. [4] for an illustra- 
tive case. 

At the critical density ratio the value of 6 is about 0.2 

(21) (cf. Fig. [2]), while the Nusselt number reaches a minimum 
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1 



Fig. 2. Determination of the thickness 6 of the stagnant 
zone, (example for Ra* = 10 7 , Le= 0.01, Pr= 0.1). The 
steep curve shows the heat flux in the stagnant zone as 
a function of the assumed value of <5, the shallow curves 
the corresponding flux in the overturning zone. Intersection 
points are possible values of S for stationary heat flow. Top: 
density ratio R p = 4. For the higher density ratio of 5.5 
(bottom), there is no value of the thickness for which the 
two match (see text). 



value of the order of a few. At this density ratio, the model 
predicts a jump from an overturning state (Nut > 1) to a 
purely diffusive state Nut = 1- An example is shown in Fig. 
|4j The presence of this jump suggests that the behavior of 
the system near the maximum density ratio needs a closer 
look than is possible with the present model. This question 
is explored with numerical simulations in ZS13. 

The maximum on R p can also be read as a minimum on 
the Rayleigh number. For an observed density ratio of 4 at 
Le= 0.01 for example, one should not expect to see double- 
diffusive layering with layer thickness less than corresponds 
to a Ra* of 10 6 (from Fig. [3j). I n terms of the reference 
length do defined below (EqT]42"|), the layers must have a 
thickness d/dg > 45, for this combination of R p and Le. 



4.2. Solute flux 

The net solut e N u ssclt num ber c an b e expressed in terms of 
Nu T ■ Using |[l8}-((^J and Q-Q, this yields the simple 
expression 



Nu s = 1 



1 



i? p Le 



1/2 



(Nu T - 1) 



(28) 



The stagnant zone of finite thickness included here thus 
leads to the same relation between Nus and Nut as the 
simpler model in S92 (cf. eq. 10 above). There is a differ- 
ence, however, in the relation between Nut and Ra* , an d 
in the presence of a maximum density ratio (cf. Figs. 3|4 1. 
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Fig. 3. Maximum value of the density ratio, as a function 
of the modified Rayleigh number of the layer, for Le = 0.01 
and Le = 0.1. 
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Fig. 4. Dependence of Nusselt number and thickness of the 
stagnant zone on density ratio R p in a range around the 



critical value, for Le = 0.01, Ra* = 10 (R 



p max 



.82). 



4.3. Asymptotic behavior 



The Nusselt number is, with ([19j), (|27j): 
Nu T = Nu T P 

where 

P 



l-Q 1 - Le 1/2 R p 
1 - Q/R P 



1 - Le 



1/2 



(29) 



(30) 



and Nut is given in terms of Ra* by (15 1. For Ra* 3> 1 



and with R p not close to R pmax , the stagnant zone is thin 



S <C 1. Ra* as a function of R p and Le is then, with (14) 
(p7): 



Ra* 



Ra* P. (Ra* > 1, R p < R p 

max J 



(31) 



The solute Nussclt number follows from (|28|), 
Nu s 



Nuq 



RpLe 1 ' 2 ' 



(Ra* > 1, R p < i? pmax ) (32) 



It is instructive to compare 5 to the thickness of the 
boundary layers Dt, D§ of the overturning zone. These 
can be written in terms of the Nussclt numbers Nut and 
Nus since at the boundaries with the stagnant zone the 
fluxes are carried by diffusion. For the thermal boundary 
layer for example this implies 



D T /(d-d s ) = l/Nu T . 



(33) 
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Using (23), (27) 



Nu T = (R P /P - 1)/S, 
hence in the limit 5 -C 1: 

D T R P 



Dt R p 
~d~ ( ~P 



1 



d 1 - 



(34) 



(35) 



Unless Q is close to unity, the width of the stagnant zone 
is thus of the same order as the boundary layers of the 
overturning zone, and the distinction between the two be- 
comes somewhat academic. At large density ratio, or under 
the fixed heat flux conditions discussed in the next Section, 
however, the distinction becomes more significant. 



5. Fixed heat flux conditions 

In laboratory situations and theoretical analyses it is usual 
to consider the temperature difference or temperature gra- 
dient as given and the heat flux or Nusselt number as the 
object to be determined. In some natural systems, however, 
the conditions under which double diffusive convection oc- 
curs are closer to ones where the heat flux is the given 
quantity. In the east African volcanic lakes, for example, 
the heat flux imposed by the influx at the bottom of the 
lake is probably more of a given than the temperature at 
the bottom of the lake. The same is the case in semicon- 
vective zones of stars. 

The temperature difference across the layers adjusts to 
an imposed heat flux. Since both the density ratio and the 
Rayleigh number depend on AT, neither of these can be 
used as control parameter any more. This requires a change 
of perspective on the problem. 

The stratifications of both temperature and solute 
change with time, under the effective transport properties 
of the semiconvective process. In the limit of low solute dif- 
fusivity, the time scale for changes in the solute profile is 
long compared with the time scale on which the tempera- 
ture profile adjusts. In this thermally quasisteady state, the 
solute profile can be taken as fixed. Assume therefore that 
the mean solute gradient dS/dz is given. As new control 
parameters, use the heat heat flux F and layer thickness d. 

To make the imposed heat flux practical as control pa- 
rameter, measure it with respect to a reference flux Fq that 
can be expressed in terms of the solute gradient. Take for 
this the diffusive heat flux that would be present in the 
linear temperature profile To(z) that is just marginally sta- 
ble against adiabatic overturning ('Ledoux'). If K is the 
thermal conductivity, the heat flux of this stratification is 



F = K 



dT 
dz 



and its density ratio is 



dS , dT 



R po = P-r-/a 



dz dz 



(36) 



(37) 



Since we have assumed that the stratification is marginally 
stable, R 



p 



1. The reference heat flux is thus 



F = K 



0dS 
a dz 



(38) 



The actual layered state has a Nusselt number Nut and a 
density ratio R p > 1; by definition of the Nusselt number, 
its heat flux is 

AT 



F = Nu T K 



With AS= ddS/dz: 



F/F = Nut^| =Nu T /i? P . 



(39) 



(40) 



To replace the Rayleigh number as control parameter, note 
that it can be written in terms of solute gradient and den- 
sity ratio as (using [T]) 



Ra* 



The quantity 



do = (g/^ = (« T /^ 



(41) 



(42) 



is a characteristic length scale of the problem: the distance 
over which temperature diffuses on the buoyancy time scale 
of the stable solute gradient, N^ 1 = (g/3 dS/dz)' 1 / 2 . The 
Rayleigh number can then be written as 



1 

R„ 



Ra* = — (d/do) 4 . 



(43) 



The Nusselt number in (40) is a fun ction of R a and Ra*, 
which can be evaluated from Eq s. ( | 15|22|23|24" ) in Sect. 
[4] above. Together, eqs. (40) and (|43p thus determine Ra* 
and R p , and other quantities of interest, as functions of the 
new control parameters d/do and F/Fq. An example of the 



dependence on d/do of Ra* 
in Fig. [| for F/F Q = 3. 



R , the Nu's, and S is shown 



5.1. Asymptotic dependences for fixed heat flux 

The limiting case d/do 3> 1 (corresponding to Ra* 3> 1), 
has a pleasingly simple form. As expected from the discus- 



sion in Sect.|4](cf. Fig.l3|, the density ratio approaches its 
maximum value Le- 1/:r in this limit. With (Eol and S 



the Nusselt numbers approach the value 

Nu s wNu T !«Le" 1/2 f/Fo {d/do > !)• 



(44) 



For Nus , this is in fact a good approximation also at lower 
values of d/do, a s Fig. [5] (lower left) shows. The correspond- 
ing effective solute diffusivity becomes 



>eS = Nu S K S « (k s kt) 1/2 F/Fq, 



(45) 



in agreement with the classical 'geometric mean of diffusiv- 
ities' scaling. The thickness of the stagnant zone is related 
to Nus by eq. (22). In the present limit, es <C 1, hence 



<5 « Nu s ~ 



Le 1/2 (F/fb 



{d/do » 1) 



(46) 



The relative thickness of the stagnant zone is thus small 
for low Lewis numbers or at high heat flux. In the oppo- 
site case of modest Le and conditions closer to marginal, it 
can be a significant fraction of the layer thickness, however. 
This may be the explanation for the relatively large thick- 
ness of the stagnant zones observed in lake Kivu (Schmid 
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Fig. 5. Dependences on layer thickness d, for fixed heat flux conditions {F/F n = 1.5). In reading order: Rayleigh number, 
density ratio, Nusselt numbers (dashed for solute) and thickness of the stagnant zone. 



et al. 2010). The heat flux measured there corresponds to 
Nusselt numbers of order 2. A stagnant zone thickness 5 
of order 30% is observed, larger than in other natural cases 
like the double-diffusive steps under the arctic ice sheet (e.g. 
Timmermans et al. 2008 and references therein). Though 
the present analysis does not apply directly to this case be- 
cause the assumption Pr < 1 does not hold (Pr ss 7 for wa- 
ter) , the thickness of the stagnant zones in lake Kivu may 
be an indication that it is actually close to the marginal 
state to be expected at imposed low heat flux conditions. 



5.1.1. Thickness of the stagnant zone 

As Fig. [5] shows, the density ra tio a pproaches its maximum 
in the limit d/d Q > 1 (cf Sect. 01): 



R p w Le _1/2 



1-Q<1. 



(47) 



The actual value of 1 — Q as a function of the parameters 
d/do and F/Fq is a somewha t co mplicated expression in- 
volving the coefficients a, b in ( 15 ) (it will not be needed in 
the following). 

Comparing S with the thermal boundary la yer thickness 
£>t of the overturning zone then yields, from ( 35 ) : 



Dt Lc 



'1/2 



- 1 



1-0 



(48) 



In this asymptotic case with fixed heat flux, the stagnant 
zone is thus much wider than the boundary layers of the 
overturning zone (but still thin compared with the layer 
thickness d). Its presence manifests its elf in the approxi- 
mate equality of Nut and Nus (eq. 44 1 as opposed to the 
original estimate (J8|. 



6. Astrophysical conditions 

For the compressible gas in a stellar interior, an incompress- 
ible approximation (Boussinesq) can still be used for flows 
on length scales small compared with the pressure scale 
height and time scales short compared with the sound travel 
time over this length, provided two factors are taken into 
account (e.g. Massaguer & Zahn 1980). First, the adiabatic 
lapse rate (dT/dz) a has to be subtracted from the temper- 
ature difference driving the flows. The modified Rayleigh 
number for a layer of thickness d is now 



Ra* 



^d i (— - — I 
kL dz dz 



(49) 



where a = \/T for the ideal gas with constant ratio of 
specific heats assumed here. In terms of the logarithmic 
gradient V = dlnT/dlnp: 



Ra* 



ft(v-v.), 



where H is the scale height of the gas pressure p, 
H = dz/dlnp. 



(50) 



(51) 



This length scale is not present in the Boussinesq case. An 
equivalent of the quantity d used in Sect. [5] still plays a 
role as well, but in the astrophysical context H is a more 
conventional choice as reference length. Ra* can then be 
written as 

Ra*=/ 2 (|) 4 (V-V a ), (52) 
where the dimensionless number /, 



/ = ( 5 f 3 /4) 1/2 , 



(53) 
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typically a large number, is the ratio of the thermal diffu- 
sion and free fall time scales over a pressure scale height. 
The dimensionless layer thickness d/H will be the control 
parameter replacing d/do above. 

A second difference concerns the radiative heat flux, 
since it is driven by the temperature gradient itself rather 
than the potential temperature^] This affects the definition 
of Nusselt number, as well as the choice of control param- 
eter specifying the heat flux. Define a thermal conduction 
constant k (a function of the opacity and the thermody- 
namic state of the gas), in terms of the the radiative heat 
flux F T : 

F T = fcV. (54) 

Following standard notation, represent the total heat flux 
F by the radiative gradient V r : 



F = kV r . 



(55) 



V r will be the astrophysical control parameter specifying 
heat flux. The ratio of density changes due to composition 
and temperature changes upon adiabatic displacement in 
the stratification, i.e. the density ratio, is 



Rp - v^v: 



(56) 



where 



V M = dln/x/dlnp, 



(57) 

and the \x the mean atomic weight per particle. Define the 
Nusselt number Nut through 



F = feV a + Nu T fc(V-V a ), 



(58) 



so that it measures the part of the heat flux that is due to 
(only) the superadiabatic part of the temperature gradient, 
instead of the total heat flux. This has to be clearly distin- 
guished from the ratio V r /V of heat flux to radiative heat 
flux in the stratification, which one might consider the more 
logical definition of a Nusselt number. The latter includes, 
however, a flux which is irrelevant for the convective flow 
(the radiative heat flux due to the adiabatic part of the 
stratification); it w ould c ont ribute even if there is no flow 
at alPj From (55), (56), (58 ) we obtain a relation between 
Nusselt number and density ratio for our astrophysical, im- 
posed heat flux conditions: 



-Rp = Nu T . 



(59) 



To complete the model, a prescription for the Nusselt num- 
ber as a function of the control parameters is needed. For 
this, we assume that the layer thickness d is small enough 
that an incompressible approxim atio n is valid, so the re- 
sults of Sect. |4] can be used. Eqs. ( 18 1— ( 26 1 then determine 



Nuj as a function of Ra* and R„. Together with R p from 



(56) and Ra* from (50), this forms a set of equations for 
the superadiabaticity V — V a as a function of the control 
parameters d/H and V r . 



2 Note that in the presence of a solute the superadiabaticity 
V — V a is not equivalent to the entropy gradient any more. 
Instead, it measures the gradient of potential temperature, the 
temperature on adiabatic displacement, in pressure equilibrium, 
relative to a given reference level. 

3 The significance of this distinction is not apparent in some 
of the numerical work on semiconvection, e.g. Biello (2001). 



6.1. Asymptotic results 

This is a somewhat implicit algebraic problem, but the lim- 
iting case equivalent to d/do S> 1 in Sect. [5] again has a 
simple form. This limit corresponds to 



f(d/H) 4 = 



gd 4 



» 1, 



(60) 
(61) 



d » z = (4#/s) 1/4 , 

i.e. d large compared to the length l on which the thermal 
diffusion time scale equals the free fall time over a scale 
height. As before, R p tends to its maximum, R p ps Le -1 / 2 
in this limit, so that 



Nun 



Le" 1 / 2 , 



while the superadiabaticity follows from (56): 

V - V a « Le 1 / 2 ^. 
The solute Nusselt number follows as in Sect. [5] 

Nus ~ Nut, 
so the effective solute diffusivity is 

Kscff = (K S K T ) 1/2 (V r - VJ/V,, 



(62) 



(63) 



(64) 



(65) 



the same as in the simpler model of S92 and ZS10 (apart 
from a factor involving radiation pressure). It is indepen- 
dent of the layering thickness d, within the range of validity 
of the approximations, 



lo < d < H. 
The relative thickness of the stagnant zones is 

<J « 1/Nus = Le 1/a V M /(V r - V B ). 



(66) 



(67) 



For Le <C V r — V a the stagnant zone is thus thin compared 
with the layer thickness (but still significantly thick er tha n 
the boundary layers of the overturning zone, cf. Sect. 5.1.1 ). 



7. Evolution of the layer thickness 

The main uncertainty of any model for semiconvection is 
the thickness d of the double diffusive layers. In a stellar 
interior the layer thic kness h as little influence on the quan- 
tities of interest (eqs. 63 65), but it can become important 
for the structure of the star when layer thickness becomes 
macroscopic, i.e. comparable with the pressure scale height. 

In geophysical and laboratory cases, it is observed that 
d evolves secularly, on a long time scale compared with the 
thermal diffusion time. Layer thickness is therefore a quan- 
tity that cannot be discussed independently of the history 
of the system. The process has not been studied very ex- 
tensively (but see Wirtz & Reddy 1979, McDougall 1981, 
Young & Rosner 2000, Ross & Lavery 2009, and the coffee 
table experiment in ZS13). The layer thickness increases by 
a process of merging of neighboring layers. Two mechanisms 
are observed: vanishing contrast, and drift of interface po- 
sition. This is illustrated schematically in Fig. [6j 

Which of these two dominates, why, and at what rate 
the merging takes place appears to be not well understood. 
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s 
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Fig. 6. Merging of two neighboring layers in a double dif- 
fusive staircase. The solute concentration can change by 
exchange between the two until the contrast disappears 
(dashed), or one can grow in thickness at the expense of 
the other (dotted). 



From the model presented here, however, an estimate for 
the merging rate which is independent of this uncertainty 
can be derived, as well as the resulting layer thickness as a 
function of time. 

Both merging mechanisms involve the redistribution of 
solute between neighboring layers. For example, if an inter- 
face between two layers transmits solute somewhat faster 
than average, the contrast between them decreases. The 
time scale r on which such changes take place is limited by 
the effective solute diffusivity, hence is of the order 



d 2 /n s 



off- 



The rate of merging is then 

dlne?/d£ m 1/r m Kg e g/d 2 , 



so that 



(2Kse ff i) 1/2 , 



(68) 



(69) 



(70) 



if Kseff is constant in time. The layer thickness is thus pre- 
dicted to increase as the square root of time. This agrees 
qualitatively with the dependence inferred empirically in 
some laboratory experiments. In most of these, however, 
the layering is induced by lateral heating of the solute gra- 
dient rather than vertical heating, hence is not directly 
comparable with the astrophysical case. Wirtz & Reddy 
(1979) for example find an initial square-root-like depen- 
dence, which saturates when the thickness reaches the sep- 
aration between the lateral boundaries. 



In terms of a merging time scale t m , (70) becomes 



d 2 /2 Ks , 



a- 



(71) 



For a quantitative check, compare this to the measurements 
from lake Kivu, where the observed layer thicknesses are 
0.3-0.5 m, the relative interface thickness S of order 0.4 
(Schmid et al. 2010). With the diffusivity of C0 2 in water, 
Kg = 2 10~ 9 m 2 /s, the effective solute diffusivity is 



Kscfi = Nu s k s ~ k s /5 « 710 9 m 2 /s. 



(72) 



For a layer of thickness 0.4 m this predicts a merging time 
of order 8 months. This agrees with the observed time scale 
for changes in the layering, of the order of several months. 
This comparison has to be taken with a grain of salt, of 
course, since the present analysis is not strictly valid for 
the Prandtl number of water. 



8. Summary 

The theory presented expands on the previous analyses in 
S92 and ZS10. It improves on these by including the effect 
of a stagnant zone of finite thickness. That is, the region 
over which heat and solute are transported by diffusion is 
not limited to just the boundary layers of the overturning 
interior of a double-diffusive step, but can in principle be an 
arbitrary fraction of the layer thickness. The need for this 
extension arose from observations of geophysical examples 
of thermohaline layering and results from numerical exper- 
iments. 

A simple 2-zone model consisting of a stagnant and an 
overturning zone, and using an experimental fitting formula 
for convective heat transport in the overturning zone pro- 
duces a clear physical picture for the dependence of the ef- 
fective transport properties of double diffusive layered con- 
vection (Sect. 4]). The results predict the existence of a max- 
imum to the density ratio for which such a layered state can 
exist. It is of the order R pmax « Le~ 1,/2 = (kt/^s) 1 ^ 2 and 
approaches this value from below with increasing Rayleigh 
number (or layer thickness). 

Of special interest is the astrophysical case where the 
heat flux rather than a temperature difference is given. In 
this case the dependence of effective solute diffusivity on 
solute stratification and heat flux has the simple form ( 65 ) . 
This is the same as before in ZS10 and S92 (except for 
the effect of radiation pressure on the equation of state not 
included here) . Due to the presence of the stagnant zone the 
value of the superadiabaticity (63) differs from that in S92, 



ZS10. For practical conditions in a stellar interior, however, 
V — V a is so small that its exact functional form makes little 
difference for the temperature stratification. The thickness 
of the stagnant zone is small relative to the layer thickness, 
as a consequence of the low value of the Lewis number. 

The main conceptual difference with respect to ZS10 
and S92 is the existence of a maximum to the density ratio. 
In the astrophysical case of imposed heat flux it plays only 
a rather implicit role, however. The differences are more 
significant when conditions are not in the astrophysically 
relevant limiting case. In numerical simulations in partic- 
ular, which are necessarily much closer to marginal condi- 
tions for double diffusive layering, the stagnant zone and 
the maximum density ratio have a substantial effect on the 
results. A comparison with such simulations is presented in 
the companion paper ZS13. 

Semiconvection is only one of the potential mixing pro- 
cesses in stars. Rotation induced mixing and magnetic pro- 
cesses are likely to be relevant as well, and could actually 
be more effective. 

Acknowledgements. The author thanks F. Zaussinger for extensive 
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